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ABSTRACT 


Experiments are carried out with various time and space 
differencing schemes applied to the barotropic primitive 
equations using both real data and a particular stream func- 
tion which is an analytic solution to the nondivergent baro- 
tropic vorticity equation. With both types of data, there 
were some significant differences in the forecasts produced 
by the various schemes. Replacement of the widely used leap 
frog (centered) scheme by others which eliminate or lessen 
some of its inherent errors at the expense of more computer 
time or storage appears to be justified at such time when 
computer capacity no longer restricts operational use of 


these more time consuming schemes. 
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fe CNDTRODUCTION 


In selecting a particular differencing scheme for appli- 
cation to numerical weather prediction, a major considera- 
tion is computer time and storage requirements. Therefore, 
the selected scheme is a compromise between verification ac- 
curacy and computer limitations. At present, the leap frog 
(centered) scheme is widely used in numerical weather predic- 
tion and will be used as a standard for purposes of compari- 
son. Present computer limitations and methods of verifica- 
tion tend to preclude seeking more accurate solutions based 
on mathematical considerations alone. The purpose of this 
study is to look beyond present computer capabilities and 
examine several other differencing schemes which reduce or 
eliminate some of the inaccuracies in the leap frog scheme. 
This scheme has three time levels and, therefore, a spurious 
computational mode arises from the finite differencing 
scheme, This problem can be eliminated by the application 
of two level time schemes such as the Euler backward or 
Runge Kutta methods which have on computational mode. A 
fourth order space differencing scheme is applied to reduce 
truncation error and improve phase speeds. Should any of 
the schemes produce a significant change in forecast rela- 
tive to the leap frog method, it must be determined if the 
difference represents a sufficient increase in accuracy to 


“warrant the added computer time and storage requirements. 





If, PROGNOSTIC EQUATIONS 


The equations of motion and the continuity equation for 


a barotropic fluid in the momentum flux form are 


9(uh) , 2 9 (uuh), 92 9 (uvh) . aes 
dt ia ox m me dy m as! eee (1) 
0 (vh) Zeon Kv uh) 290 (vvh) Cars dh 

wm «UCT ™ é(Cbx a Sy a ee mee | ee (2) 
ah , 2,8 (uh), 3 (vh) 

e) ie ae ee m Oy m : (3) 


where h is the depth of the fluid (or height of a pressure 
surface in the atmosphere), g is the force of gravity, f is 
the coriolis parameter and u and v are horizontal velocity 


components. 


A. SPACE DIFFERENCING 
1. Second Order Scheme 
Vth ievemeand HM piven at each grid point, the 


differential equations may be approximated by 


2 el ae) j mee 
—m OUT ON ee + ae ) (4) 
vhs 4 j ili 
= Aan 3 ene 
ee oe amet 


Cmeecul). . 
dt *J 


vy hie) he 

us eaeste: 

eee cS Y m 
vh, fal eee 
ne ee ? me 


+ (u ay 


= roy en 


- mgh..(h 


eee a 4/24 = Fy, 


13 
where j increases with x and i decreases with y. (See Fig. 


1) Similar expressions are used for the other two equa- 
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tions. This method of space differencing, devised by A. Ara- 
kawa (see, for example, [1]), prevents nonlinear instability 
in the continuous time case (sometimes referred to as the 
semi-discrete equations). 
2. Fourth Order Scheme 
Fourth order space differencing which takes into 
account a greater number of grid points and reduces trunca- 


tion error approximates uh tendencies by 


2 uh 


peecub), . . 1 m __it2, 3 
at by] OS Lf Re are ag Sn se) 
uh, j uh, _» j ou 
cara cS ae aae ? woe 
en VWs oy ae MD le 
vh vh 
2 age 
Ga a ae Ba =" ) 
vh, Foe are 
ee Oh) am =P ia )] + Soo ONDE 
STehiG SES - he» 3 ) jaa S 
where F . is the second order space difference. 
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Doe tiME DIFFPERENCING 
ese Pwesetee suns Vil, OF Ws t 1S time; n is time step 
number and F is -: in the following time differencing 


schemes: 


ieee ap oh s0)e 


gers 2 y= = aie 


2. Runge Kutta 


AU At F[u", ndAt] 


AU, = At F[U"+ AU,/2, (nt5) At] 


lie 





“AU 


_ n Z die 
AU, = At F[U + oar e (n + 5) Ave | 
AU, = At Pe Se Ug 5 Ge Ne 
n+1 


n 1 
U + glAu, + 2AUq + 2AU, + AU, J 


Cc 
tt 


ou Euler Backward 


* 
teow = CA UF 


* 
Uo -ue C= SCOAtSOéSF 


* 
Denotes trial values. 
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Pine DOUNDARY CONDITIONS 


pee EARTICULAR STREAM FUNCTION 
The boundary conditions in this case are applied to an 
initial stream function which is an analytic solution to the 
nondivergent barotropic vorticity equation and after each 
meme Step during the forecast. | 
Meeevorcn— Olen boundary Conditions 
For the 63x63 grid a wall exists between y=l1 and 


y=2 and also between y=63 and y=62 (See Fig. 1). No flux 


is permitted across this rigid boundary. This is acconm- 
plished by setting Lea e Be, and eg a thereby 
making the average Ye 7 QO along the wall. The other para- 
b 
meters are set equal across the wall. Summarizing, the 
boundary conditions are 
North . South 

63 ia Vx ,62 el aa Yx,2 

u = 

eo. = “x ,62 “x1 7 rx 32 

ce 62 Oeil 7 ame nD, 
where @ is the geopotential. 

2. East-West Boundary Conditions 


The East-West boundary conditions are cyclic as 
suggested by the function plotted in Fig. l. Summarizing, 


the boundary conditions are 


3 





East West 
A = A A = A 
ey. O25 Ory. eS 


where A represents parameters v 2 and @ 
X,Y X,Y x,y 
3. Variation for Fourth Order Space Difference 
In order to apply fourth order space differencing 
and still maintain the boundary conditions of the original 
prid, it was convenient to create an additional outside row 


of points and relabel as a 65x65 grid. (See Fig. 2) The 


boundary conditions are summarized as follows: 


North South 
vx,65 VS 62 Yx,1 7 Vx 44 
A ge Cretan is ; 
Vv =-Vv V = -Vv si and @ edad 
x,64 x./03 > ees xa xy. 
65 Seetay: Saal 7 ae 
ees) x 63 7, Panes 
East West 
A = A A = A 
ee B25 y 65,y 4,y ee represents ete 
A = A A = A Vv and ; 
2,y 63,y 64,y 3,5y X sy ° Py 


B, ACTUAL 500 MILLIBAR DATA 

When the equations were integrated using real data, the 
"spongy" boundary conditions of the FNWC model were applied. 
After each time step, the prognostic results were modified 


by setting 


uh = (1-k')uh + k' uh, 
vh = (1-k')vh + k' vh, 
ea eeleio Vinee kk’ he 
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where the ,. subScript indicates initial values. Below lat- 
itude four degrees north set k'=4. Above 17 degrees north 
set k'=0. Between four degrees north and 17 degrees north, 
Bere sis linear variation in k' from zero to one. This re- 
stores the variables to their initial values south of four 
degrees north, while north of 1/7 degrees north, full varia- 
bility is permitted. There is no restoration north of 17 


degrees north. Between, there is partial restoration. 


15 
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Stream function. Ihe dimensionless sym- 


bols x and y are used for notation con- 
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PMMictiidttimom tne —OoxOS Expansion of the grid in 
Pomeloimmoraer sto apply fourth order space dif- 
ferencing over the same area. The dimensionless 
symbols x and y are used for notation convenience. 


Figure 2 
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A. A PARTICULAR STREAM FUNCTION 

MeooxOossduare grid as in Fig. 1 and as modified in 
Fig. 2 for fourth order space differencing is used to 
minimize the number of changes required to adapt programs 
to the FNWC 63x63 hemispheric grid. The boundaries in 
latitude extend from 30 to 90 degrees. Grid length is 
determined by dividing the geographic distance represented 


by this latitude difference into 62 equal parts. 


B. FNWC 63x63 HEMISPHERIC 

In Conjunction with the actual 500 millibar data, the 
FNWC 63x63 hemispheric grid is used. This is a square re- 
peesentation Of a hemispheric polar stereographic grid in 


which the grid distance d is determined by 


t 2 O 
dt _ 40/2 a sin 60 
m eo San ¥ 


ii 


d ) 


: ; : : O 
where m is the map factor, d' is the grid distance at 60 


latitude (381 kilometers), and sin y is defined by 


PD ie 
sin y = 2f3 tt = where Ro = (I-I yee (J-J ie 
973.71 + R u : 
Plo -wacpDlitary @rid point. 
Sane PaerceepiOot mes tOtrthie north pole = 32,32. 
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A. A 


voor LNL ULAR SCONDITIONS 


PARTICULAR STREAM FUNCTION 


LMmoGEecrt ee tOmobpkain 4a Smooth Initial stream function 


E1eld, 


WOrticity equatdioneis 


where cis 
ection and 


HO, Old (CL 


Here, Yo 


~~ 
~~, 


v 


O 


where v 
O 


x 


is determined 


an analytic solution to the nondivergent barotropic 


used, namely 
sine (y - gy cos cai (x-ct) 
i Z L 
- x 


the phase speed and the wavelength in the x dir- 


the half wavelength in the y direction are equal 


Ly -—Soldj. Lhe inttial stream function is 
sin —— (y - oS cos ei x 
i. 2 L 
va. = 


by using the relation 


Veale 
Oo x 
25 


is the initial specified amplitude of the meri- 


dional wind speed derived from 


V= 


Using the initial field thus obtained, 


equation 


oy 
Ox 


27 x] 
= : 
x 


= 2u 


: TT d 
L Vv, [sin 7 (y- 5) cos 
x y 


the linear balance 


Woo = Vm 4s We oan 


which in finite difference form may be 


V*¢..= 


fad 
+) 1j 


2 
V Vig 


written 


itt 
sl 


ill 
+ lvees Nee Pes Vea ij 


Meg) 





where 
et 2-157 6 a, ger* O99, g-1- 80a, 


Pe y= CS: 


The preceding equation may now be solved as a Poisson 
equation for the geopotential field. Thus, 
2 = 
V ° 5 ae 
The initial guess for the relaxation procedure is 


5 = f vo where £ is a mean coriolis parameter. 


Subsequent "guesses" are made according to the relation 


(nt1). g(a) 4g (a) 
ij 


p 


where 


a) ne Vat) 
soe = a(V Oe = hee 


and @ is a relaxation coefficient. 

teoone h = o/¢, computations for the initial height field 
are completed. Additionally, a constant height is added 
£O each grid point value in order to eliminate negative 
heights, which would otherwise cause computational insta- 


bility. (See Fig. 3) 


B. ACTUAL 500 MILLIBAR DATA 

The initial height field is obtained from the FNWC 500 
millibar analysis. The linear balance equatian is solved 
for the stream function as a Poisson equation using W= o/£ 


Somat irst fuess. 


ZO 





Differencing RMSE SboOr Position Mister or Central 


Scheme (Nautical Miles) Height (Meters) 
Leap Frog 97 30 

Leap Frog Fourth 76 31 

Order Space 

Runge Kutta pe 31 

Euler Backward 88 25 


Table of root mean square errors in position and central 
height of the three low centers at 24, 48, and 72 hours 

as forecast by various differencing schemes (on actual data) 
in comparison to the corresponding analysis. 


TABLE 1 


Latitude 





Latitude 


Initial height field as produced by a parti- 
cular analytic stream function on a 63x63 
oridweeOlly= =the central contours of the three 
features are illustrated. 


Figure 3 
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Ve “RESULTS 


The results shown in the following figures were obtained 
from the computer subroutine mapping of the 63x63 height 
field obtained by the various differencing schemes using 10 
Minute time steps. One grid squares are used in the figures 
Simply as a convenience to provide a ready reference scale 


to compare positions of low centers. 


A. PARTICULAR STREAM FUNCTION 

Results for the forecasts made using a particular 
Stream function for the initial values are summarized in 
Fig. 4 for central height and relative position. The fore- 
cast results for the individual differencing schemes are 


Preccembed in greater detail in Fig. 8 through Fig. 12. 


B. ACTUAL 500 MILLIBAR DATA 

The results of forecasts made on three low systems from 
the 500 millibar FNWC analysis are summarized in Fig. 5, 
Fig. 6 and Fig. 7. Although only low systems are illustra- 
ted in this study, similar results were obtained from the 
movement and intensification of highs. The time frame of 
this study begins with the 1200Z March 24, 1966 500 milli- 
bar heights and terminates with the 1200Z March 27, 1966 
forecast. Low number one (See Fig. 5) was located over 
Somieneoweoi bettas Low number two (See Fig. 6) and low nunm- 


ber three (See Fig. 7) were located over Southwestern 


eae 





Canada and Labrador respectively. Root mean 
computations for position and central height 
by each differencing scheme when compared to 


mreeanalysis are contained in Table l. 


twee COMPUTER REQUIREMENTS 
Although peak program efficiency was not 


poe the application of the various numerical 


Square error 
as forecast 


the correspond- 


attempted dur- 


schemes, con- 


sistency in application was carefully maintained. Using 


the leap frog scheme as a standard, the leap 


PROC we Our en 


order space differencing and Euler backward schemes require 


appmominately a jo percent increase in run time. The Runge 


Kutta scheme requires approximately a 170 percent increase 


in run time and a 25 percent increase in computer storage. 
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ere Oe ——— ve ene 


Central Height (Meters) 


Forecast Dh rene Gi n\re ey aa ae 
x Leap Frog Sie 4 383 383 
oLeap Frog Fourth 

Order Space 383 384 383 
aRunge Kutta | Boe 385 380 
vEuler Backward 333 384 Sieh 
@Initial Central Height --- 383 Meters 


Relative position and central height of a 
particular stream function feature as fore- 
cast by various differencing schemes. 


Figure 4 
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an 24 hr 
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| | 206 Nautical Miles 
eal 
One. grid square 


EN teen ef 


Differencing Scheme Central Height (Meters) 


Forecast ZAC Nr 7-2 NY. 
xLeap Frog 4930 4940 4990 
reap Frog Fourth 

Order Space 4940 4960 5030 
vyRunge Kutta 4930 4940 4990 
vEuler Backward 4940 4940 4990 


Corresponding 
FNWC Analysis ©5020 ©5090 @5110 


@iInitial Height 
FNWC Analysis 5040 


Forecast and analyzed relative positions and 
central heights for an actual data 500 milli- 
bar feature designated low number one. 


Figure 5 
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Corresponding 
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@Initial Height 
FNWC Analysis 5100 


Forecast and analyzed relative positions and 
Ceomenbatemerehtis for an actual data 500 milli- 


bar feature designated low number two. 


Figure 6 
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Forecast and analyzed relative positions and 
ConeralewiereittsceLOr an actual data 500 milli- 
bar feature designated low number three. 


Figure 7 
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Vite CONCLUSIONS 


The various differencing schemes were applied to the 
smooth analytic solution in order that differences in fore- 
casts produced by these schemes would be more readily dis- 
tinguishable than with actual data. However, the choice 
as to resolution was probably too high in this study (61 
grid lengths per wavelength) and shorter wavelengths are 
Meteto be tested. Differences produced in location and 
central values were not conclusive when compared to the 
results produced using actual data. A reduction in resolu- 
Buoweror the analytic initial field (approximately six grid 
lengths per savellena@ny would possibly produce a more mean- 
Pork separation in the accuracy produced by the various 
schemes. 

The negligible increase in accuracy (four percent based 
On a root mean square error comparison) of the Runge Kutta 
scheme over the leap frog scheme does not justify the ex- 
tremely large increase in computer time and moderate in- 

a accetmmeconpubker Storage. Both the leap frog fourth 
order space and Euler backward differencing schemes show 
promise for improved accuracy with a moderate increase in 
computer run time (approximately 35 percent). The leap 
frog fourth order space differencing scheme gave a 24 per- 
cent increase in position accuracy. The Euler backward 


TIwememprodieced sd Nane percent increase in position accur- 
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acy and a 17 percent increase in central height. With both 
the Euler backward and the leap frog fourth order space 
schemes, the 35 percent increase in computer time could be 
reduced somewhat by more efficient programs. The small 
statistical sample of pressure systems examined in this 
study does not permit firm conclusions to be drawn regarding 
the various schemes tested. However, there is reason to 
believe that the fourth order space Ae eee eta of the 
nonlinear advection terms will reduce the phase error in the 
prediction of pressure systems and this method should be 
tested further. Williams [2] has shown that the pressure 
force need only be differenced with second order accuracy 
which permits a somewhat larger time step without computa- 
tional instability than when fourth order differences are 


used throughout. 
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One feature in an initial height field pro- 
duced from a particular stream function 
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Figure 8 
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Forecast height in meters of central values and 
positions of centers and 390 meter contours for 
epaeewctlar stream function by the leap frog 


scheme. 
Figure 9 
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Forecast height in meters of central values and 
positions of centers and 390 meter contours for 
a particular stream function by the Runge Kutta 


scheme. 
Figure 10 
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positions of centers and 390 meter contours for 


a particular stream function by the Euler back- 
ward scheme. 
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One grid square 


Forecast height in meters of central values and 
positions of centers and 390 meter contours for 
Gepant@eular Stream function by the leap frog 
POtreimotder. Space differencing. 

Figure 12 
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